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Abstract 



J2 
o 

We extend quantum kinetic theory to deal with a strongly Bose-condensed 
atomic vapor in a trap. The method assumes that the majority of the vapor is 
not condensed, and acts as a bath of heat and atoms for the condensate. The 
condensate is described by the particle number conserving Bogoliubov method 
developed by one of the authors. We derive equations which describe the 
fluctuations of particle number and phase, and the growth of the Bose-Einstein 
condensate. The equilibrium state of the condensate is a mixture of states 
with different numbers of particles and quasiparticles. It is not a quantum 
superposition of states with different numbers of particles — nevertheless, the 
stationary state exhibits the property of off-diagonal long range order, to the 
extent that this concept makes sense in a tightly trapped condensate. 
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I. INTRODUCTION 

In two previous papers on quantum kinetic theory, which we shall call QKI [1] and QKII 
[2] in this paper, we focussed on the theory of a Bose gas in which there was 

a) No more than a small amount of Bose condensate, so that the interactions between 
the particles were not sufficiently strong to produce a significant modification to the 
excitation spectrum of the Bose gas; 

b) No trapping potential, so that the system was spatially homogeneous. 

Such a treatment gives insight into the the statistical aspects of quantum kinetics, but 
cannot be applied to the condensates presently being produced [3-8] , in which there is very 
strong condensation of an alkali metal vapor in a rather tight trap. 

In this work the focus will be on strongly condensed trapped systems, with a methodology 
based on that of QKI, but somewhat simplified in order to exhibit the major features of 
the interaction of a Bose-condensate with a vapor of non-condensed particles. The basis of 
our method is to divide the condensate spectrum into two bands; the condensate band, in 
which the presence of the condensate significantly modifies the excitation spectrum, and the 
non- condensate band, whose energy levels are sufficiently high for the interaction with the 
condensate to be negligible. 

In Fig. 1(a) we depict the kind of trap for which this division is very clear — a rather wide 
trap, in the center of which is a small tight trap, with a few energy levels. In contrast, 
in Fig.l(b) we depict a harmonic trap, in which the division is not so obvious. In fact, 
however, the three dimensional nature of the trap, in which the density of states increases 
quadratically with energy, also yields a relatively small number of isolated levels at low 
energies, which rapidly approaches a continuum, so the dramatic difference between the 
two traps is more apparent than real. Current experiments all use harmonic traps, but the 
"ideal" trap of Fig. 1(b) would in principle be better adapted to studies of the condensate 
itself. 

The dynamics of the non-condensate band is thus well described by the kind of formalism 
used in QKI, with modifications to take account of the fact that we have also introduced a 
trapping potential Vy(x). On the other hand, the dynamics of the condensate band is well 
described by the energy spectrum which arises from Bogoliubov's method [9], as adapted to 
apply to the case of particles in a trap [10]. This kind of dynamics has been well studied by 
many workers in the past two years [10,11,22], and the equation of motion for the condensate 
wavefunction is now generally verified both theoretically and experimentally to be the time- 
dependent Gross-Pitaveskii equation [12]. 

The aim of this paper is to develop and use methods to describe the interaction between 
the condensate band and the non-condensate band. The results of our investigations can be 
summarized as follows: 

i) The non- condensate band is taken to be fully thermalized, in the sense that at any 
point u, k in phase space there is a local particle phase-space density F(u, k), and the 
space time correlation functions of this phase space density have a weakly damped 
classical form. 



ii) The dynamics of the condensate band, in interaction with the thermalized non- 
condensate band, is given by a master equation (50a-50f), in which the interaction 
terms are of three kinds: 

(a) Two particles in the non-condensate band collide, and after the collision, one of 
these particle remains within the condensate band; and there is also of course the 
corresponding time-reversed process. 

(b) A particle in the non-condensate band collides with a particle in the condensate 
band, and exchanges energy with it, but does not enter the condensate band. 

(c) A particle in the non-condensate band collides with one in the condensate band, 
and after the collision both particles are left in the condensate band; and again 
there is also of course the corresponding time-reversed process. 

hi) By using the particle-number-conserving Bogoliubov method [18] developed by one of 
us, the expressions which arise in the master equation can be re-written in a form in 
which the processes are described in terms of 

(a) Creation and destruction of particles; thus if N is the total number of particles 
in the condensate band, we introduce operators A\ A which change N by ±1. 

(b) Creation and destruction of quasiparticles, which are to be regarded as quantized 
oscillations of the whole body of N particles in the condensate band. 

Thus there is a clear distinction made between creation and destruction of particles 
(N — > N ± 1), and the creation and destruction of quasiparticles, which are mere 
excitations, and not particles. This distinction follows from the particle number con- 
serving Bogoliubov method [18], and has not been made clear in nearly all 1 earlier 
treatments of Bogoliubov's method, which blur this distinction, and give the mislead- 
ing impression that the Bogoliubov excitation spectrum can only be obtained if strict 
conservation of particle number is abandoned. We emphasize here that the particle 
number conserving Bogoliubov method gives exactly the same results for the excita- 
tion spectrum as the usual version; however, because particle number conservation 
is maintained, it is possible to specify both wavefunctions and energy eigenvalues for 
states with a definite energy and a definite number of particles. 

iv) It is possible to split off the largest terms in the master equation, and approximate 
the description by a single basic master equation for N, the number of particles in the 
condensate band, and thus arrive at a description of the growth of a condensate from 
the vapor. The very simplest version of this description, in which all fluctuations are 
neglected, takes the form of a single differential equation 

N = 2W + {N) {(l - e [M-M(JV)]/fe B T^ N+1 ^ (1) 



1 A notable exception is the work of Girardeau and Arnowitt [19] 



in which W + (N) is a certain kinetic coefficient, for which we can give an approximate 
form, and \x and fi(N) are respectively the vapor chemical potential and the chemical 
potential of the condensate wavefunction for N particles, as given by solving the time- 
independent Gross-Pitaevskii equation. 

v) Detailed descriptions of the dynamics of fluctuations can be given. The stationary 
solution for the condensate density operator is a mixture (not a quantum- mechanical 
superposition) of states with different numbers N of particles. For each such state, 
the wavefunction of the condensate is £jv(x), the solution of the the time-independent 
Gross-Pitaevskii equation for N particles. 

iv) A description of both phase diffusion of the condensate and of the decay of excitations 
can also be given, and these will appear in a subsequent paper. [13]. 

The description given in this paper is a simplified one, in which the non-condensate band 
is treated as stationary, so that it does constitute a genuine "heat-bath" , characterized by a 
fixed temperature T and a fixed chemical potential \i. This is not an essential feature of the 
basic quantum kinetic description, which, in its full form, treats the two bands as having 
their own dynamics as well as being coupled together. The full quantum kinetic theory will 
be presented in a subsequent paper, QKV. It turns out to be a rather intricate process to 
derive such a description, even though the essential result, as used in this paper, is that the 
non-condensate band is described by an equation of the Uehling-Uhlenbeck form [15], with 
correlation functions of a thermal form. 



II. DESCRIPTION OF THE SYSTEM 

A. Hamiltonian and notation 

As in QKI, we consider a set of Spin-0 Bose particles, described by the Hamiltonian 

H = H kin + H I + H T , (2) 



in which 



and 



tf kin = |rf 3 x^t( x ) (-^V 2 ) ^(x) (3) 

Hj = - /rf 3 x / > rf 3 xV t (x)^ t (x>(x-x')^(x / )^(x). 

(4) 

The potential function w(x — x') is a c-number and is as usual not the true interatomic 
potential, but rather a short range potential — approximately a delta function — which repro- 
duces the correct scattering length. This enables the Born approximation to be applied, and 
thus simplifies the mathematics considerably. [17] 

The term Ht arises from a trapping potential, and is written as 



H T = J rf 3 x\/ T (x)^( x )^(x). (5) 

Thus, this is the standard second quantized theory of an interacting Bose gas of particles 
with mass m. 



B. Condensate and non-condensate bands 

When there is strong condensation in a trap, the lowest energy level of the trap develops 
a macroscopic occupation, and modifies the other energy levels in the trap. If the trap is 
quite tight, the lower energy levels are reasonably well separated, as is indeed experimentally 
the case [3-5]. The concept of a "tight" trap is however not absolute; a trap may be tight 
for the lower energy levels, but quite broad for the higher energy levels. The essential 
distinction is the density of states, which increases like E 2 for a harmonic trap, leading to 
very closely spaced levels for quite moderate energy. In addition, one must also note that 
the modification of the energy levels induced by the condensate affects only the lower levels 
significantly. 

The procedure which follows from these is to divide the particle states into a condensate 
band Re and non condensate band Rnc, and to perform a corresponding resolution of the 
field operator in the form ^(x) = 0(x) + ipNcfe)- The separation between the two bands 
is in terms of energy; we choose to define the division in such a way that the energies of 
all particles in Rnc are sufficiently large that there is essentially no effect on them by the 
presence of the condensate in the lowest energy levels. We will describe Re in terms of 
discrete trap energy levels, while Rnc will be taken as being essentially thermalized, with 
correlation function given by the usual equilibrium formulae. 

1. Separation of condensate and non- condensate parts of the full Hamiltonian 

We can now write the full Hamiltonian in a form which separates the three components; 
namely, those which act within Rnc only, those which act within Re only, and those which 
cause transfers of energy or population between R c and Rnc- Thus we write 

H = H NC + H + H ItC (6) 

in which H N c is the part of H which depends only on ipNC, H is the part which depends 
only on 0, and that part of the interaction Hamiltonian which involves both condensate 
band and non-condensate band operators is called Hjc, and can be written 

H hC ee ff« + Hi% + H®, (7) 

where the individual terms in Hj c are the terms involving operators from both bands, which 
cause transfer of energy and/or particles between Rq and Rnc- We call the parts involving 
one 4> operator 



H {1) - 
n i,c — 



rf 3 X / rf 3 x'M(x - x')^]vc( X )V ; Arc( x/ )V ; A r c(x)0(x') 

+ /d 3 x /rf 3 x'M(x-x / )0 t (x)^j vc (x / )^c(x)^c(x') 



The parts involving two operators are called 

rr(2) __ 



| rf 3 x | rf 3 x' M (x - x / )^j vc (x)0t(x')Vivc(x)0(x / ) 
+ J rf 3 x I rf 3 x' M (x - xO^x^^xOV^x^x') 

+^ I rf 3 x I rf 3 x' M (x - x')^ c (x)^ c (x')0(x)0(x / ) 

+i J rf 3 x I rf 3 x' M (x - x')0 t (x)0t( x > J vc(x)^c(x') 



The parts involving three operators are called 



#(3) = 
n I,C — 



rf 3 X / rf 3 x'M( X - X / )0 t (x)0 t (x , )0(x)^ JVC (x / ) 
/Vx /"d 3 x'u(x - X')^] vc (x)0 t (x / )0(X)0(X / ). 



(8) 



(9) 



(10) 



III. DERIVATION OF THE MASTER EQUATION 

Using the separation into condensate band and non-condensate band operators we can 
now proceed to develop a master equation. We have developed a full quantum kinetic 
description, in which the atoms in Rnc are treated by a wavelet method as in QKI, and the 
atoms in Re in terms of the wavefunctions corresponding to the energy levels of the trap, 
modified by the presence of a condensate in the lowest level. This leads to a description 
which is very comprehensive, but one whose complexity would obscure the essential aspects of 
condensate dynamics which are the subject of this paper. This description will be published 
separately. 

The present description will assume that Rnc is fully thermalized, and thus the atoms 
in these levels will be treated simply as a heat bath and source of atoms for the levels in Re, 
which provide a description of the dynamics of the condensate treated as an open system. 
The simplest assumption is that p^c has a thermodynamic equilibrium form given by 

PNC=ZW{ ^ ) (11) 



This gives rise to relatively simple equations for the kinetics of the condensate under the 
restriction that the density operator for the non-condensed atoms is clamped in this form. 
A more practical thermalization requirement is that of local thermalization. This is defined 
by requiring the field operator correlation functions to have a thermal form locally, and to 
have factorization properties like those which pertain in equilibrium. Time-dependence of 
Pnc can be included as well, as long as this happens on a much slower time scale than than 
those which arise from the master equation which we derive. The precise nature of these 
local equilibrium requirements will be specified in Sect.IIIA5 

A. Formal derivation of the master equation 

The derivation of the master equation follows a rather standard methodology. We project 
out the dependence on the non-condensate band by defining the condensate density operator 

as 

p c = Tr NC (p) (12) 

and a projector V by 

Vp = p NC ® Tr NC (p) (13) 

= Pnc ® Pc, (14) 

and we also use the notation 

Q=l-V (15) 

The equation of motion for the full density operator p is the von Neumann equation 

p = -^[H NC + H + H I<c ,p] (16) 

= (C NC + £ + Ci, c )p. (17) 

We use the Laplace transform notation for any function f(t) 

/(*)= / e- st f(t)dt. (18) 

Jo 

We then use standard methods to write the master equation for v(t) = Vpit) as 

sv(s) - v(0) = C v(s) + VC IiC v(s) 
+VjC i , c [s - C NC -Co- QCi^QC^cvis). (19) 

In this form the master equation is basically exact. We shall make the approximation that 
the kernel of the second part, the [ ] _1 term, can be approximated by keeping only the 
terms which describe the basic Hamiltonians within Re or Rnc, namely the terms Cnc and 
Cq. We then invert the Laplace transform, and make a Markov approximation to get 

v(t) = (£o + VCi,c)v(t) 

+ \vC ItC f Q dr exp {(£ + Ci,c)t} QC I>c v(t)} . 

(20) 

There are now a number of different terms to consider. 



1. Forward scattering terms 
These arise from the term VjCi y cv(t), and lead to a Hamiltonian term of the form 

-^forward = / (fix. / rfVufx - x') 



forward 

xTr JVC |^] vc (x)^ c .(x)0t(x / )0(x / ) 
+^c(x)^c(x')0 t (x / )0(x)l (21) 



The terms involving the operators represent essentially two effects. The first is the effect 
of the mean density of the condensate band on the energy levels of the non-condensate band. 
By definition of the non-condensate band, this effect is negligible. They also represent an 
effect of the average non-condensate density on the condensate. This is also small, but can 
be included explicitly in our formulation of the master equation. 

In (21) one can see clearly the "Hartree" term in the second line, and the "Fock" term in 
the third line. However in the remainder of this paper we shall use the delta function form 
of the interaction potential w(x — x') — > uS(x — x'), and the difference between these terms 
and the second and third lines of (9) disappears. None of our results will depend heavily on 
the validity of this approximation, but the formulae do become considerably more compact. 

2. Interaction between Re and Rnc 

We now examine the term involving H IC as defined in (8), which contain one 0(x) or 
</>'(x); explicitly, the term in Hamiltonian can be written 

Hfy = f d 3 x Z 3 (x)0 f (x) + h.c. (22) 

In this equation we have defined a notation 

Z 3 (x) = wV4c( x )V^c(x)^Arc(x). (23) 

Substituting into the master equation (20), terms arise which are of the form 

-^ J rf 3 x J d 3 x' J™ dr Tr {z 3 (x)Z 3 t ( X / , -r)p NC } 

x0t( x )0(x',-r) Pc7 (t). (24) 

Here we have used the notation 

Z 3 (x, t) = e iHNct/h Z 3 (x)e- iHNct/h (25) 

0(x, t) = e iHot/h (j)(x)e- iHot/n (26) 

We can now transform the term involving the 0's into something more tractable, though 
perhaps more abstract, by expanding in eigenoperators of the commutator [H , ]: Thus we 
write 



0(x) = ]TX m (x) (27) 

m 

where the operators X m (x) are eigenoperators of H 

[H ,X m (x)] = -he m X m (x). (28) 

Here e m may be positive or negative, but m is a label which uniquely specifies the eigenvalue 
e m . The quantities e m give the possible frequencies corresponding to transitions between 
condensate energy levels, and will be called the transition frequencies. 

3. The nature of the operators X m 

To compute the transition operators X m , one needs to know all of the eigenstates of Hq. 
Thus, if these eigenstates are called \l,c), with eigenvalue Ei, then 



X m (x) = ^5(e m ,^-^)|/',c')(/',c , |0(x)|/,c)(/ 



c . 



(29) 



This expression is of a formal nature, since the full calculation of the energy eigenstates 
and the matrix elements of the field operator is a formidable task. However, by using the 
modified Bogoliubov method of [18], this can be made into a tractable task, as will be 
explained in Sect. Ill C. 

4- Random phase and rotating wave approximations 



Because the two operators <p and (jy both turn up, there will be a double summation 
^ mm tX\ n X m i. To get a master equation it is necessary to ensure that only the m = m! 
terms exist. This is normally achieved in quantum optical situations via the rotating wave 
approximation, which effectively eliminates the unwanted m ^ m! terms on the grounds 
that they are rapidly oscillating, unlike the m — w! terms. Perhaps more compelling in 
the situation here is the fact that the diagonal (non-rotating) terms have a consistent sign, 
while the off-diagonal terms — of whiche there are very many — do not, and will thus tend 
to cancel. This is effectively a random phase approximation. Hence we will keep only the 
diagonal m = m! terms in what follows. 

5. Correlation functions of Rnc 

The terms involving Z^,Z\ are averaged over pnc, which is quantum Gaussian, but is 
modified by the presence of the interaction. A full treatment of these correlation functions 
requires the quantum kinetic theory of the non-condensate band, and this will be done in [14]. 
However this treatment is quite intricate, and yields a result which, in the approximation 
that the non-condensate is regarded as a thermal undepleted bath, is really quite standard 
and well-accepted. This takes the form 

9 



i) We may make the replacement 

<Z 3 (x)4(x\-r)> 

x(^ nc (x)^ nc (x>,-t)) (30) 

(Terms involving (ip* NC ('x.)ipNc(' x -') etc. do arise in principle, but give no contribution 
to the final result because of energy conservation considerations.) 

ii) The time-dependence is needed only for small r, and in this case we can make an 
appropriate replacements in terms of the one-particle Wigner function F(u, K) 



J (» + ^)l>Nc(u-^, 



(27f) 3 Jr nc 



f d 3 KF(u,K)e- iK - v e luj(K ' u)T 



(31) 



|)^c(u-|, 



'/'JVC ( u + - ) V'jvc ( u - -> ~ T 



^ / d 3 K [F(u, K) + 1] e *K-v e -MK,u)r 



(27f) 3 JRatc 



(32) 



where 



hu(K,u) = ^- + V T (n). (33) 

Since the range of all the K integrals is restricted to Rnc, it is implicit that in all 
integrals that hu(K, u) > Er. 

iii) This approximation is valid in the situation where F(K, u) is a smooth function of its 
arguments, and can be regarded as a local equilibrium assumption for particles moving 
in a potential which is comparatively slowly varying in space. 

6. Final form of the term 
Using these factorizations, the term (24) becomes 

- /rf 3 xrf 3 x'^(-)(x,x / ,e m )^Xt i (x)X m (x')pc;(t) (34) 

in which 
0("> (u + ^, u - |, w) = —^2 / d 3 K, J d 3 K 2 J d 3 K 3 [F(Kx, u) + 1] [F(K 2 , u) + 1] F(K 3 , u) 

xe i(K 1+ K 2 -K 3)l P (wi + ^ - ^ - w) (35) 

= G(") (u + |, u - |, w) - iG, H (u + |, u - |, w ) (36) 

10 



In the above we use the notation 

SP(x) = - f" e~ lXT dT (37) 

TV JO 

= 8(x) - --. (38) 

■K X 

The imaginary part in (36) gives rise to level shifts, which we shall neglect in what follows, 
both for simplicity and because they are likely to be very small. 



7. Terms involving two 4> operators 

Of these terms, only those involving on and one $ can yield a resonant term, cor- 
responding to scattering of a non-condensate particle by the condensate. The effect of the 
non- resonant terms would be a shifting of the energy levels of Rnc, an d by definition of the 
non-condensate band, this is negligible. 

Thus we arrive at an approximation to H\ c in the form of a term 

HJfy « 2 j d 3 x Z2(x)0 f (x)0(x) + h.c. (39) 

where 

Z 2 (x)=^tvc(x)V4cW- (40) 

The term analogous to (24) in the master equation becomes 



4 poo 

tf Jo 



2 , dT0t( x )0(x)0(x',-r)0t(x', 



xJd 3 xJd 3 x f (Z 2 (x)Zi(x.', -r))p c (t). (41) 

We can write now 

0(x)0t( x ) = J2 (]TX m (x)Xt(x)5(e m - e n ,p r )) (42) 

r \m,n / 

= E^(x). (43) 

r 

Where f3 r represents the frequency associated with U r (x). Inserting these into (41), we get 

~¥ r dT I d ' X I rfV < Z 2(x)Z 2 t (x', -r)> 

xEe^(x)f/ r (x> c (t). (44) 



11 



8. Terms involving three <f> operators 
The part of the interaction Hamiltonian can be written 

Hi% = J rf 3 x^ c (x)0t(x)0t( x )0( x ) + h .c. (45) 

The master equation term this time takes the form 

x dr0 t (x , ,-r)^(x',-r)0(x , ,-T) / 9 C (t). (46) 

Jo 

We can write now 

0t( x )0(x)0(x) 

= E f E ^(x)X n (x)X p (x)5(e m - e n - e p ,#)) 

r \m,n,p / 

(47) 
^EK(x). (48) 



Inserting these into (46), we get 



"h? L dT I rf3x / rf3x ' ^ JVC ( x ^c( x '> ~ T ^ 



X 
r 



J2e^ T VXx)Vr(x')pc(t) (49) 



B. Final form of the master equation 

By carrying out the procedures used to derive the term (34), we can finally get the master 
equation: 



PC it) = - 



l - J d 3 x 0t (x) (-^- + V T (x) + 2uJ d 3 K F(K, x) + ^t( x )0( x )^ ( x ) , pc 



(50a) 



/ d 3 x f dV f ^ G (+) (u, v, e m ) (2X m (x)p c Xt l (x') - p c X r t „(x')X m (x) - Xl(*)X m (x)p c ) (50b) 

+ 5]G(-)(u,v,e m ) (2Xt,(x)p c X m (x') - p c X m (y!)Xl(y) - X m (x')^(x) Pc ) (50c) 

m 

+ ^M(u,v,/3 r ) (2C/ r (x)p c (7 r t (x') - Pc C/ r t(x')C/ r (x) - l/t(x / )E; r (x)pc) (50d) 

r 

+ £)£<+> (u,v,#) (2T/ r (x)p c V;t(x') - pc^x')^) - ^(x'^x)^) (50e) 

+ £)£;<-)(u,v,#) (2F r t(x)p c y r (x') - p c T/ r (x')V;t(x) - V r (x')y r t(x)pc) J • (50f) 



12 



In this equation we use the notation 

u=(x + x')/2 (51) 

v = x-x', (52) 

and define the quantities E^\ M, G^ by 

G<+> (u, v, to) = " /" d 3 Kx /" d 3 K 2 y d 3 K 3 e iAK ^ v F(K u u)F(K 2 , u) [F(K 3 , u) + l]5(Aoo 123 - w) 

(53a) 
G(-) (u, v, w) = —^—2 y d 3 Kx y d 3 K 2 y d'W^^-^Ki, u) + l] [F(K 2) u) + 1]F(K 3 , u)6(Acj 123 - w) 

(53b) 
M(u, v, W ) - -^p y d3K i y (PKie'^-K'^FQd, u)[F(K 2 , u) + l]<5(wi - c^ 2 - «) (53c) 

fl(+)(u,v,u;) = "' 2fi2 Jd 3 K 1 e iK ^F(K 1 ,u)6(u; 1 -w) (53d) 

£H(u,v '" ) = 2(^¥/' l3KielKrv|f(Kl ' u) + 11% ^ ) - (53e) 

Here we have used the notation 

Ldi = Lu(Ki,u) (54) 

AK 123 = Ki + K 2 - K 3 (55) 

Aa;i23 =uj ± +u 2 - uj 3 . (56) 



1. Relation between backward and forward rates 

If we consider the case that pnc corresponds to the thermodynamic equilibrium form 
(11), we can set 

F(K, u) -»• N(K, u) (57) 

= {exp [(hu(K, u) - n) /k B T] - l}' 1 (58) 

and it follows, for example, that 

F(Ki,u)F(K 2 ,u)(F(K 3 ,u) + i) e -(M-^ 123 /fc B T) 

= (F Ql (u) + l)(F Q2 (u) + l)F Q3 (u) (59) 

so that in this case 

G (+) (u, v, w) = e (M-MA B T G (-) ( U) V) ^ ^ 

M(u, v, u) = e- hu}/kBT M{u, v, -w) (61) 

#' (U, V, w) = e b-*»)/k B T E {-) ( u> v? w ) (g 2 ) 
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C. Use of the Bogoliubov approximation for the condensate 

The master equation in the form (50a-50f) gives an accurate treatment of the internal 
dynamics of the condensate band, with an approximate treatment of the coupling to the non- 
condensate band. However in order to use it in practice, it would be necessary to compute 
the full spectrum and eigenstates of the condensate, and for this it is necessary to have some 
approximate treatment of the condensate. The Bogoliubov approximation is the natural tool, 
but in its conventional form there is the difficulty that the Bogoliubov eigenfunctions do not 
have a well defined number of particles. In [18] a treatment of the condensate spectrum was 
given which combined together the Gross-Pitaevskii equation and the Bogoliubov spectrum 
to give approximate eigenfunctions with an exact number of particles. We shall use this 
treatment in our paper, but one should bear in mind that all that is needed is a description 
of the energy levels for each number N of atoms in the condensate band. In particular, 
the work of Girardeau and Arnowitt [19], based on the pair model, also presents such a 
description in the case of an untrapped condensates, which could be ganeralized to this 
situation [20], and this would probably yield a more accurate description. 

This treatment can be summarized as follows for the situation here, where we wish to use 
it for the condensate band, which will involve in principle some adaptation to take account of 
the fact that </>(x) is expanded in basis functions which belong to Re- However the changes 
involved are so minor that we will not consider them in detail here. The major effect is 
caused by the fact that the commutator of the condensate band field operators is not a delta 
function, but takes the form 

[0(x),0t( x ')] = 4( X)X ') (63) 

where g(x, x') is the restriction of the delta function to the condensate band. This is in- 
dependent of the amount of condensate, since the boundary between the condensate and 
non-condensate bands is defined in terms of the unperturbed eigenstates. 
The operator 0(x) is written in the form 

0(x) = a 6v (x) + J2 Ofc6( x ) (64) 

k 

where £at(x) is the ground state wavefunction for the condensate in a situation where there 
are N particles in the condensate band. All quantities on the right hand side are therefore 
implicitly functions of N, and this expansion is thus different for each N. The wavefunctions 
£fc(x) together with £jv(x) form a complete orthonormal set. 

The essence of the method is first to find an approximation to 0(x) which is valid when 
0(x) acts on a highly condensed state, in which most of the particles are in the state 
represented by the wavefunction £jv(x). Such a highly condensed state is written |no,n), 
where n = {n^} represents the vector of occupation numbers of the non-condensed modes 
within the condensate band. 

If iV = hq + J2k n k, the state can be rewritten in a form which uses the total number of 
particles N rather than no and we will call this \N, n). One now defines the operators 

A\N,n) = \/N\N -l,n) (65) 

M = A^A (66) 

oik = —rj= r <Aa k (67) 
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If N is very large, we can write approximations valid in the space of fixed N 

ao*A(l- ^) (68) 

w A. (69) 

The simple form (69) is adequate in almost all situations, except the actual computation of 
the eigenf unctions. The inverse of (67) is also approximately given by 

Aak tnn\ 

a k ps . (70) 

VNTT 

and the field operator is then approximated by 

' T E^(xK) (71) 

1 k ) 

X(x)l • (72) 




(The corresponding formulae in [18] use y/N, rather than \/N + 1, which is what appears 
most naturally in the computations. The difference is of order 1/N, which is better than the 
accuracy of the expansion, but the \/N + 1 form gives much more elegant formulae because 
of the cancellation with other terms — hence its use here.) Finally we can show that the 
phonon operators a k , and the operator A behave approximately like a set of independent 
annihilation operators; i.e., 

[a k ,a\] ps S k ,k'- (73) 

It is then found that if 6v( x ) represents the ground state wavefunction, the following are 
true. 



1. Gross-Pitaevskii equation 

The ground state wavefunction £jv(x) must be a solution of the time-independent Gross- 
Pitaevskii equation 

-^V 2 6v(x) + VHx)6v(x) + AH6v(x)| 2 £v(x) 

= /i W 6v(x). (74) 

In fact, because of the modified commutator (63), the precise equation is given by the 
projection of the right hand side of (74) onto the condensate band being equal to the left 
hand side, but there will be very little practical difference introduced by this modification. 
It is, roughly speaking, equivalent to solving the partial differential equation on a grid of 
spacing similar to the lowest wavelength in the non-condensate band. This is in fact what 
is done in any practical procedure — the only modification is that we now specify that the 
grid not be finer than required by the restriction to the condensate band. 
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2. Approximate Hamiltonian 

The procedure can be put in the form of an expansion in inverse powers of iV 1 / 2 by 
formalizing the requirement that the number of particles N be large and the interaction 
strength u be small. To do this one sets 

u = u/N (75) 

and then develops the approximation procedure as an asymptotic expansion at fixed u in 
inverse powers of N 1 ' 2 . Using this procedure the condensate band Hamiltonian can be then 
approximated by 

H = NHi + H 3 (76) 

where 

Wi = -^/A^(x)V 2 Mx) + |d 3 xG(x)VHx)6v(x) + ^/d 3 x|6v(x)| 4 , (77) 

= E (N)/N (78) 

which is a c-number, and 

W 3 = -|^yd 3 xx t (x)V 2 x(x)+yd 3 xx t (x)T/(x)x(x) 

+ J d 3 xj|(6v(x) X t(x)) 2 + ^(G(x)x(x)) 2 +x t (x)x(x)(2S|^(x)| 2 - m )|-^|d 3 y|^(y)| 4 . (79) 

3. Chemical potential 

The condensate chemical potential fJ,(N) is determined by the requirement that 

rf 3 x|^(x)| 2 = l (80) 

and this will give a different fi(N) for each N: It can be shown that 

fi(N)^E (N+l)-E (N). (81) 

4- Bogoliubov transformation 
The Hamiltonian 7i 3 can be diagonalized by a Bogoliubov transformation of the form 

a k=Yl C kmb m + Yl S kmb\n ( 82 ) 

m m 

and here b m is a quasiparticle destruction operator. We can then write 

X (x) = ]T (^(x)6 m + <4(x)^) (83) 
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with 

Pm (x) = E4 m 6(x) (84) 

k 

fc(x) = Ei(k(4 ( 85 ) 

k 

The diagonalized Hamiltonian is then written 

H 3 = huj g (N) + J2 he m (N)blb m . (86) 

m 

The practical task of determining the values of the necessary quantities is not trivial, but 
it is at least perfectly well defined. It does not matter, in principle, what functions £fc( x ) 
are chosen as long as they are a complete ortho normal set in the space orthogonal to £at(x). 
However, a choice of appropriate ^(x) may make the practical problem of diagonalizing 
H3 somewhat more straightforward. There have been a number of papers [21-23] in which 
numerical methods of diagonalizing 7i 3 have been presented 

5. Relationship between ground states for N and N + 1 particles. 
The operator A does depend on N because a is defined by (65) to be 

a o (AO = y>xG(x,iV)0(x), (87) 

where the explicit dependence on N of both £^ and ao has now been written. 

This means that a Q (N)\N,n = 0) is a state with N + 1 particles in the wavefunction 
corresponding to ground state of the N particle state; it is not the ground state for N + 1 
particles. We will introduce the notation \M)n to mean a state with M atoms in the 
wavefunction corresponding to the ./V particle ground state. 

It can be shown that the operator B^N) which connects the iV and N + 1 ground states 
through 



B*(N)\N) N = VNTT\N + 1) N+1 (88) 

B(N-1)\N) n = Vn\N-1) n _ 1 (89) 



is given approximately by (to order 1/N) 



B\N)^(l + ^=J2rlal\aUN). (90) 



where 



j3_ d£jVV 



r k = Njd 3 x^^Ux,N)- (91) 

This means that the field operator expansion (72) now takes the form 
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0(x)«£(JV)(&v(x) 

+ 7FTT^ (x) "^ (x) ? rfcafc 9- 

(92) 

It is particularly interesting to see that the corrections to the quasiparticle term are of the 
same order of magnitude as the original terms; the correction is thus quite significant. 

D. Transformation of the master equation 

The Bogoliubov approximation for </>(x) (72), together with the expression of x( x ) as 
in (83) has the advantage that it gives, almost directly, the operators X m (x) which were 
implicitly defined in (27). 

1. Transition operators in the master equation 

We can see that there are three kinds of transition operators as defined in (27), but their 
form depends on the number of condensate band particles N in the state on which they act. 
Thus, we can write the action of on an arbitrary state in the N basis as 

4>(x)\A) N = (x (x) +]T{X m (x) + X+(x)}) \A) N 

(93) 
in which we have the correspondence, with eigenvalues and eigenf unctions, as follows: 

Operator Representation Eigenvalue 

X (x) B(N - l)6v-i(x) _ -/i(iV-l) (94) 

X-(x) B(N-l)b m f m (N-l,x)/y/N -%e m {N - 1) - fi(N - 1) (95) 

X+(x) B(N-l)blg m (N-l,x)/y/N he m {N - 1) - n{N - 1) (96) 

X (x)t B(NYC N (x) fi(N) (97) 

X m (x)t B(N)W m fc(N,x)/y/W+l he m (N)+fi(N) (98) 

X+(x)t B(Nyb m g* m (N,x)/VNTT -%e m {N) + fi(N) (99) 



in which 



f m (N, x) = p m {N, x) - ^iv(x) ]T r kC km (100) 

k 

= Pm(N,x) (101) 

-iV^(x)/rf 3 x'p m (iV,x')^^, (102) 

g m (N,x) = q m (N,x) - £ N (x) ^r k s km (103) 

k 

= q m (N,x) (104) 

-iV^(x) J dV q m (N, x')^^, (105) 



The importance of the corrections can now be seen — they are of the same order of magni- 
tude as the direct quasiparticle terms, and represent the generation of quasiparticles as a 
consequence of the changing shape of the condensate. 

As given above, the choice of how to write the operator 0(x) depends on the value of 
N, and in a master equation we have to consider the action of several <fi operators on the 
left or right of a density operator, which may contain off-diagonal terms. To see what is 
involved, suppose we take a term like <ftp<f>, and consider a term in p like \Ni)n 1 {N 2 \n 2 - It 
is soon clear that the use of the N± expansions on the right and the N 2 expansion on the 
left leads to an evolution operator which is extremely complicated. However, in practice we 
will not be interested in density operators which are far from diagonal. This means that we 
will instead write off-diagonal terms in the form \N — u)n(N + v\n, and for such terms the 
expansion appropriate to N will be the correct one. An off diagonal term of this kind is an 
eigenfunction of the two-sided operator N, defined by 

Np = l[N,p] + (106) 

= l[B*(N-l)B(N-l),p] + (107) 

An arbitrary density operator can then be expanded in the form 

PN = Y. r v\ N -»)N(N + v\N- (108) 

We now treat the problem of expanding master equation terms of the kind 

Term 1 = 20(x)p JV t (x') - [p N , t (x')0(x)] + , (109) 

Term 2 = 20t( x )p w 0(x') - \p N , 0r(x')0 t (x)] + . (110) 

We use the the expansion of the field operator as follows; we work only to lowest order, that 
is we omit all quasiparticle terms, for which the working is analogous: 

0V)0(xW - t (x')^-i(x) J B(iV - l)p N 

- &-i(xO£t(JV - l)6r-i(x)£(JV - l)p N . (Ill) 

The use of the B^(N — 1) follows because B(N — 1) takes a term like \N — u)n — > |iV — 1 — 
v)n-i, and then B^(N — 1) takes it back to \N — v)n Using this methodology, we then find 
that 

Term 1 = 6v-i(x)G_i(x) [lB(N - \) Pn B\N - 1) - [p N , B\N - 1)B(N - 1)] + } , (112) 

Term 2 = 6v(x)G(x') {2B(NYp N B(N) - [p N , B(N)B^N)]+} . (113) 

In order to get a form in which the resolution of the density operator into terms p^ is not 
necessary we use the operator iV to rewrite (112) in a form valid for any p 

Terml = ^B {U^Mp} B ^ ~ [{&-i(*)&-i(*>} . B ^ B } + } , ( 114 ) 

Term2 = ^ {^(x)^(xOp} B - [{&(*)&(*>} , BBt ] + } ( ll $) 

In this form it must be understood that 

19 



• p is expanded as a sum of terms pn- 

• The operators B, B* are defined by 

B\A) N = B(N-1)\A) N (116) 

(A\ N B = (A\ N B(N) (117) 

B*\A) N = B\N)\A) N (118) 

(A^St = (A\ N B^(N - 1) (119) 

The terms involving quasiparticles can be similarly treated. 

2. Other transition operators 
The operators U r (x) and V r (x) are similarly able to expressed using these expressions. 

3. Master equation in terms of the Bogoliubov states 

We can now write out the master equation in terms of six transition probabilities defined 
in terms of functions R ± as 

W + (N) = R + (£ N ,p N /h) (120) 

W-(N + l) = RT(Z N ,fi N /h) (121) 

W+ + {N) = R+(f m , ( C ™ + p N )/h) (122) 

W~-(N + 1) = R-(f m , (e% + p N )/h) (123) 

W+~(N) = R + (g m , (-eft + ^)/n) (124) 

W~+(N + 1) = i?"( 5m , (-e£ + p N )/h) (125) 

The functions R ± (y,ui) are defined by 

i? + (y, w) = " 2 / d 3 x /" d 3 K 1 rf 3 K 2 d 3 K 3 d 3 k^Aa )123 (x) - o;)5(K! + K 2 - K 3 - k)FiF 2 (l + F 3 )W„(x, k) (126) 

R~(y,u;)= " /d 3 x /d 3 K 1 rf 3 K 2 d 3 K 3 d 3 k5(Au; 12 3(x)-a ) ) ( 5(K 1 +K 2 -K 3 -k)(l + F 1 )(l + F 2 ) J F 1 3 >V 1/ (x,k) 
(27r) D ?i J J 

(127) 
Here we use the notation 

w » (x - k >=(^/ d M x+ iM x -i) e ' kv (l28) 

to represent the Wigner function corresponding to the wavefunction y(x). 

We now use the explicit expression (53a) for G^ + ^(u, v,cj), and substitute for the X 
operators to get 
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50b^2B{w-(N)p}B^ - [b*B, {w~(N)p} 






blB^BL 



w- 



\N)p \ 



N 



UtMi™) 



A' 



> J + 



The reversed process leads to the term given by 

50c^ 2B ] {W + (N) P } B - [BB\ {w + (N)p} 



. 2Bt6 ,jii»} ftmB 



fc t BBt6 ,„,{l™^} J+ 



(129) 



(130) 



^. More complex terms 

The terms involving two and three condensate field operators would become very com- 
plicated unless an abbreviated notation is devised. We introduce a notation 



B = S 
Bbl = S... 



(131) 
(132) 
(133) 



so that while m is an index whose range is conventionally the positive integers, there is now 
an index 1 = 0, ±1, ±2, . . . which enumerates the full range of operators. 
We also introduce the notations 



and for the frequencies 



L o (x,iV)=£(x,A0 
L m (x,A0 = / m (x,AQ/v /iVTT 
L_ m (x, N) = (? m (x, N)/y/N + l 
W/(iV,x,k) = W L/ (x) 



p(N)/h = Q (N) 
p(N)/h + e m (N) = n m (N) 
p(N)/h-e m (N) = n_ rn (N) 



(134) 
(135) 
(136) 
(137) 



(138) 
(139) 

(140) 



In terms of these abbreviations, we can write all the two terms in terms of a set of rate 
functions 

Au 2 



Ru(N) 



I d 3 x f d 3 K! f d 3 K 2 I d 3 k f d 3 k'5(K! - K 2 - k + k') 



(27r)^ 2 7 ~~j --v w " 2 y w ~y »-^«i «. -.-, (141) 

xF(Kx, x)(l + F(K 2 , x))Wj (iV, x, k)Wj(JV, x, k')5(Acu 12 (x) - fi 7 (iV) + Slj(N)) 
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and this leads to the transformation 

50 rf - Y, hSiSj {Ru(N)p} SjSJ - [SjSJSjSl {R u (N)p}]} (142) 

The three terms require forward and backward rate functions 

R+ JK (N) = U 5 2 /* d 3 x j d 3 K! / d 3 k / d 3 k' j d 3 k"5(K! + k - k' - k") 

xF(K 1 ,x)W/(iV, x, k)Wj(iV, x, k')Wx(iV, u, k")<5 ( Wl (x) + n 7 (JV) - Qj(N) - Q K (N)) (143) 

RJ JK (N + 1) = — ^-2 / d 3 x j d 3 Ki / d 3 k j d 3 k' / d 3 k"5(Ki + k k' k") 

x(l + F(Ki,x))W/(Ar,x,k)>Vj(Ar,x,k , )WK(Ar,u,k")<5(wi(x) + f7 / (iV)-f7,7(iV)-f7K(iV)) (144) 
and then take the form 

50e - 2 (25/5J4 {i? / + JX (A r )p} S^jS] - [^SjSJ^J^, {i?^(iV)p}] } (145) 

UK ^ ' 

The reversed three <fi terms give 
50/ -> 2 {25}5j5a; {iW^)p} S ] K S]S I - [sk^SfSjS*, {iW^)p}L| (146) 

UK L ' J 



IV. PRACTICAL APPLICATION OF THE MASTER EQUATION 

The master equation using all the transition operators is a rather complex object, and it 
is wise to attempt a simplification which would give only the most significant contributions in 
order to get some insight into the basic structure predicted. The most obvious simplification 
is to take into account only the terms which are most significant when there is a large amount 
of condensate. This means that we will consider only the terms which are of highest order 
in N, and these are easily identifiable as: 

a) Those terms arising from the first lines of (129,130); 

b) The terms involving only R 00 (N) in (142); 

c) The terms involving only R 0Q0 (N) are certainly almost zero because they would involve 
the absorption of an atom from the non-condensate band into the condensate level, 
which cannot conserve energy. Thus we can probably ignore these terms completely. 

A. A master equation for the condensate mode alone 

For a first investigation we will consider only the terms from the first lines of (129,130), 
since the Rqq(N) terms cannot change the number of particles in the condensate. This leaves 
a simplified master equation for the condensate mode alone in the form 
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-J:\KP] 



+25 1 " [w + {N)p] B - [BB\ [w + {N)p] 
+2B {W'(N)p} fi f - [B ] B, {w-(N)p} 



(147) 



Here H' Q = H + 2M/d 3 K/<i 3 xi r '(K, x)0(x)t0(x) includes the mean-field correction to the 
trapping potential, which is probably very small. 

The ground state wavefunction is in practice very localized compared to the size of the 
cloud of atoms in the non-condensate band. This means that F(K,x) is very broad in 
x compared to W(£, x, k), which is comparatively sharply peaked at x = 0, so that for 
example in the integral in (129) over x we can simply use F(K, 0) and Acji 23 (0). With 
these approximations, the transition rate functions (126, 127) are 



W + (N) 



W~{N + 1) 



u 



fd 3 K± fd 3 K 2 fd 3 K 3 /rf 3 k5(K] 



K, 



K 3 - k) 



(2ufh 2 

x5(Auj U3 (0) - p(N)/h)F(K 1 ,0)F(K 2 ,0)(l + F(K 3 ,0))\£ N (k)\ 2 



'14$ 



u 



f rf 3 Ki f d 3 K 2 I d 3 K 3 f d 3 k <J(Ki + K 2 - K 3 - k) 



{2ufh 2 

x5(Acu 123 (0) - p{N)/h){l + F(Kx, 0))(1 + F(K 2 , 0))F(K 3 , 0)|| JV (k)| 2 (149) 



in which <^v(k) is the momentum space condensate wavefunction 

1 



6v(k) 



(2tt) 3 / 2 



cfxe 4k ' x 6v(x) 



(150) 



B. Thermal bath of non-condensate atoms 



In this case we assume that we can write 



3 (ftw(K,x)-/i)/fc B T 



F(K,x) ; 
from which one easily obtains (choosing Vr(0) = 0) 

W + {N) = e { ^^ N)}lkBT W'(N+l). 



fl51^ 



(152) 



i. Maxwell-Boltzmann approximation 

The integrals can be approximately evaluated by making a Maxwell-Boltzmann approx- 
imation F(K,x) ps e -(?iw(K,x)-/i)/fe s T_ Although this may be in some cases a drastic approx- 
imation, it has in its favor: 

(i): It should always give the correct order of magnitude; 
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(ii): The distribution is needed only for energies in the non-condensate band, and will not 
become infinite there as long as Er > \x\ 

(iii): It will certainly be valid for sufficiently large K. 

Using u = Aiiah 2 /m, where a is the s-wave scattering length, we get 

k B T)}' 

(153) 



nh I k B T 



Here K\(z) is a modified Bessel function — the full details of the approximations involved in 
the evaluation of W + (N) are given in Appendix B. 

2. A fuller quantum treatment 

Evaluation of W + (N) beyond the Maxwell-Boltzmann approximation involves a more 
serious consideration of the cut at Er between the condensate and non-condensate bands. 
When one makes the Maxwell-Boltzmann approximation the dependence on the "cut" at 
Er is not very strong; however, when the full Bose-Einstein form is used this dependence 
becomes more severe because of the Bose enhancement of the distribution function at lower 
energies. The correct treatment of this problem must take into account the particular 
situation being treated. For example, if we are considering the growth of the condensate, 
the condensate and non-condensate bands will not be in equilibrium with each other, and the 
dependence on Er obtained may be very misleading since it arises from the lower energies 
at which one would expect some modification of the pure Bose-Einstein distribution. One is 
led to the conclusion that this issue can only be resolved by taking into account more details 
of the dynamics of both the lower levels of the non-condensate band and of the quasiparticle 
levels in the condensate bands, which should provide the correct interpolation between the 
lower end of the non-condensate band and the upper levels of the condensate band. 

On the other hand, when we are considering the dynamics of the condensate in equilib- 
rium with the non-condensate — i.e. the equilibrium fluctuation dynamics of the condensate, 
the non-condensate distribution function must be the Bose-Einstein function, subject only 
to the levels being describable as the unperturbed trap levels. 

C. Stochastic master equation for diagonal matrix elements 

We can easily take the diagonal matrix elements of (147) and derive the simple birth- 
death stochastic master equation for these 

p(N) = 2NW + (N - l)p(N - 1) 

+2(N + l)W~{N + l)p(N + 1) 

-(2(7V + 1)W + (N) + 2NW~(N))p(N). (154) 

This is of a very standard form from which we can derive the following. 
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1. Stationary solution 
This is given by 

and using (152), this can be written 

P S (N) = exp ( J2 \M ~ »(M)]/k B T) (156) 

\M=1 ) 

<xex P ([fiN-E (N)]/k B T), (157) 

and is the result expected from statistical mechanics. We can make a Gaussian approxima- 
tion to P S {N) by expanding about its maximum, which occurs at the value N determined 
by \x — [a(N), and the variance of this Gaussian is given by 

var(iV) = k B T/^B. (158) 

The Thomas-Fermi approximation to the condensate chemical potential is fi(N) oc iV 2 / 5 , 
and when this is valid it can be seen that 

var(iV) oc N 3/5 (159) 

which can be sub-Poissonian for sufficiently large N, although this may not be achieved in 
practical cases. 

2. Order •parameter and off'- diagonal long range order 
The correlation function of relevance is 

(0t( x , £)0(x', t)) 8 = £ iVP s (iV)G(x)6v(x') (160) 

N 

» iVG(x)&v(x') (161) 

whenever var(A^) is very small compared to N, which is clearly the case for large enough 
N. However it is perfectly clear that (0^(x, t)) = (0(x',t)) s = 0, so that off-diagonal long 
range order is achieved without breaking the phase symmetry of the Hamiltonian 0(x) — > 
exp(i#)0(x). 

Of course the order achieved here is clearly not "long range", since the ground state 
wavefunction £jv(x) is only non-zero over a very small distance interval, unless the trap is 
very broad, which would no longer be within the range of validity of the present treatment, 
which assumes distinctly spaced energy levels. 
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3. Fokker-Planck equation 

For sufficiently large N we can make a Kramers-Moyal expansion [24,25] to the stochastic 
equation (154), yielding the Fokker-Planck equation 

p(N,t) = 2^| [NW-(N) -(N+ 1)W + (N)} p(N,t) 



d 
+ 



NW + (N)p(N,t)]\. (162) 



ON 
This is equivalent to the stochastic differential equation (in the Ito form) [24,25] 

dN = 2[(N + l)W + (N) - NW~(N)] dt 



+2^/NW+(N)dW(t). (163) 

Assuming once more small fluctuations, we can linearize about the mean stationary value 
N by writing N = N + n, and use the relationship between the forward and backward rates 
(152) to obtain the linear stochastic differential equation 

dn=-kndt + \ / DdW(t) (164) 

in which 

k = 2 m H N)^ a -§l (165, 

D = ANW + (N). (166) 

The characteristic relaxation time of the number fluctuations is thus given by 1/k for small 
fluctuations about the equilibrium value. 

4- The growth equation 

The rate equation arising from (154,163) for the mean number is obtained by neglecting 
all fluctuations and takes the form 

N = 2[(N + 1)W + (N)-NW-(N)}, (167) 

and when the thermal form for relation between the forward and backward rates (152) is 
taken, this becomes 

A> = 2W + (N) { (l - e (M(A0-M)A s ^ TV + l} . (168) 

This equation combines aspects of laser theory and chemical thermodynamics into one simple 
statement. The picture is analogous to that of two chemical species, "condensate" with 
chemical potential fi(N), and "vapor" with chemical potential //, which come to equilibrium 
when the number of atoms in the condensate iV satisfies the requirement 

N = \l ( e M#)-ri/*BT _ i) . (1 69 ) 
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Since it is assumed that N >> 1, this means that p(N) = p, to order 1/N. This mean 
corresponds to the modal value, iV mod , of P S (N) determined by W + (N mod — 1) = W~(N mod ) 
which from (152) requires p(N mot j) = p. 

We will call equation (168) the growth equation; it can easily be integrated numerically, 
and the results have already been published in [26]. The characteristic behavior is best 
illustrated when we start with the initial condition with N = 0, and assume that the vapor 
has sufficient atoms for us to neglect its depletion as atoms enter the condensate, so that p 
is assumed to be constant. Under these conditions the initial growth of N is given by the 
"spontaneous" term in the growth equation (168), i.e., N = 2W + (N). This gives a relatively 
slow growth, until N becomes large enough for the "gain" term proportional to N become 
significant — this term is positive if p > /i(0), and yields a dramatic net gain, giving a rapid 
rise in N which saturates as N approaches the stationary value N s , for which p = p(N s ). 
Rubidium 



D. Phase diffusion 

This section gives the simplest possible discussion of phase diffusion, for which we will 
consider the single-mode master equation in the form (147). The additional terms (142) will 
have a significant effect on phase diffusion, but for simplicity are not included here. A full 
treatment of phase diffusion will be given in [13]. 

The phase decay is given by the equilibrium time-correlation function G(x, t, x', t') = 
(0^(x, t)0(x', t')), and we can evaluate this using standard techniques [27] 

G(x,t,x',t') =C(x)£WTr{Ai{e L WAp a }} (170) 

where L is the master equation evolution operator defined by writing (147) in the form 
p = Lp. Notice that the operator expansion in terms of A is used since this preserves the 
basis state expansion; thus A' p s can be written straightforwardly as an operator of the type 
(108). 

If we define 

p'{r) = e LT A Ps (171) 

then the object we need to compute is 

Tr{Atp'(r)} = J2VN^1{N- l\p'(r)\N) (172) 

N 

subject to the initial condition 

Tr {AV(0)} = ]>>(iV|p s |iV> = (N) s (173) 

N 

Let us first compute Qn(j) = (N — l\p'(r)\N) subject to the initial condition ^(0) = 
v~N(N\p s \N). The master equation yields an equation of motion for ^(r) in the form 
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9n — — r — 9n 

+2^/N(N -l)W + (N - l)g N -i - (27V + l)W + (N)g N 

+2^/N(N + l)W~(N + l)g N+1 - (27V - l)W + (N)g N . 

(174) 

We now assume that g^ is only significantly different from zero when TV is large and carry 
out an expansion similar to the Kramers-Moyal expansion [24,25]. We first set 



2 V /TV(TV±1)^2TV±1- — (175) 

and then make the Kramers-Moyal expansion to get 

(ifi(N) W + (N) + W-(N)\ 
9N * {— 47V ) 9N 

+ 2 dj^{i NW+ ( N )9N}- (176) 

The Fokker-Plank operator in the second two lines of this equation is exactly the same as 
that if (162), but the two terms on the left of the top line cause a major difference in the 
behavior of the solutions. The simplest way to elucidate the behavior of the solutions is to 
linearize the equations in much the same way as led to the approximate stochastic differential 
equations (163-166). This means that we write approximately 

//(TV + n) « MTV) + n^p- 

= n(N) + hen, (177) 

W + (N) + W~(N) _ W + {N) + W~(N) 

4TV W 4TV 

= k. (178) 

It is necessary to go to higher order in the first equation, (177) because this term is imaginary, 
and may cause significant cancellation, which cannot happen for the second term, which is 
real. The linearized equation then takes the form, written now in terms of the function 
f(n, r), defined by g{n + TV, r) = exp[(i/i(TV)/fr — K,)T]f(n, r) of the variable n 

f(n) = ienf(n) + |- \knf\n) + ~f(n)\ . (179) 

This equation can be solved by the same method as used in [24] Sect. 3. 8. 4. The initial 
condition follows from that for g(N), and, to the degree of accuracy being used in this 
linearized treatment, can be taken as f(n, 0) = NP S (N + n). This leads to the solution for 
the characteristic function 
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f(s,T) = jdne ms f(n,r) (180) 

/ Ds 2 De i , x De 2 \ , s 

-^{-ar-wi 1 -* )-w T ) (181) 

From the initial condition (173) it follows that 

Tt{aV(t)} = e^^-^N f dnf{n,r) (182) 

= e {i ^ N)/n ~ K)T Nf{0,r) (183) 



eM^-^JVexp 



2 



2A; 3 



1 - kr - e~ kT 



(184) 



This equation shows the expected oscillatory behavior determined by the chemical potential 
(a(N), and three characteristic damping time constants. 

(i): The time constant l/«, which is the analogue of phase diffusion in a laser, and has its 
characteristic form of the noise divided by the number of particles. We call the corresponding 
time constant 

47V nh 2 N 

Tl = W+(N) + W-(N) ~ 2ma 2 {k B T) 2 (185) 

(ii): The behavior for fct < 1 of the final exponential factor is of the Gaussian form, 
characteristic of inhomogeneous broadening, exp [— (r/r 2 ) 2 /2], where 



r-2 



h 2 jd^N))- 1 

^i^iv-j (186) 



(iii): For kt ^> 1, the behavior of this term becomes exponential with a time constant 

2k 2 8ma 2 N , m . 

Ve z Tin 

This time constant is independent of temperature, unlike the other two. 



E. Stochastic master equation including quasiparticles 

Let us now take into account the quasiparticle terms which arise from the terms (129,130) 
and develop a stochastic master equation for the occupation probabilities. We can define an 
occupation probability 

p(N,n) = (N,n\p\N,n) (188) 

where n = {n m }, the set of all quasiparticle occupation numbers. The master equation terms 
(129,130) then give rise to a stochastic master equation for these occupation probabilities 
p(N, n), in the form 
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p(N, n) = 2NW + (N - l)p(N -l,n)-2(N + l)W + (N)p(N, n) 
+ 2(N + 1)W~(N + l)p(N + 1, n) - 2NW~(N)p(N, n) 
+ ]T{2n m W++(iV - l)p(N - 1, n - e m ) - 2(n m + l)W+ + (N)p(N, n)} 

m 

+ ]T{2(n m + 1)W~(N + l)p{N + 1, n + e m ) - 2n m ^"-(iV)p(iV, n)} 
+ ]T{2(n m + l)W+-(N - l)p(N - 1, n + e m ) - 2n ro W+-(JV)p(N, n)} 
+ ^{2n m W" + (iV + l)p(7V + 1, n - e m ) - 2(n m + l)W~ + (N)p(N, n)} (189) 

m 

Here e m — {. . . 0, 0, 1, 0, 0, . . .} has its only nonzero value at the position corresponding to 
the index m. 

There is no contribution to this equation from the terms (142), since these do not transfer 
particles between the vapor and the condensate, and we assume at this stage that the 30 
terms in (145,146) are negligible. 

1. Scattering of quasiparticles 

The stochastic master equation (189) is clearly incomplete, since it lacks the terms which 
would arise from the terms (142), which represent the scattering of quasiparticles by atoms 
in the noncondensate band, thereby effecting a transfer of population between quasiparticle 
levels. The inclusion of these terms will not affect the stationary solutions, but could cause 
significant changes to time dependent quantities, such as condensate growth or time correla- 
tion functions. Nevertheless we shall not pursue these aspects now, but will deal with them 
in a later publication. 

One should also note that there can also be scattering of quasiparticles by each other, 
but these terms only appear if one takes the expansion in inverse powers of yN of the 
Hamiltonian H to the next two higher orders after those considered in (76). This aspect 
will also be left to a later publication 

2. Stationary solution 
Using similar reason to that in (IV B) we can show that 

W+ + {N + 1) = eb i -rt N) - em]/kBT W—(N) (190) 

W+~(N + 1) = e [ ^^ N)+em]/kBT W~ + (N), (191) 

and from these it is straightforward to show that the stationary solution of the stochastic 
master equation (189) takes the expected form 

P S (N, n) oc exp i — S . (192) 
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3. Deterministic equations of motion 

Assume factorization of all correlations, and neglecting 1 compared to N, the determin- 
istic equations of motion, analogous to the growth equation (168), take the following form. 
First define 

n+ = 2W+ + {N) ((1 - e (M(Af) -^ +£m)/fesT )n m + l) (193) 

n- = 2W~ + (N) ((1 - e^ N)+ » +e ™ )/kBT )n m + l) (194) 



Then 



1T"m 


n m + n m 




N : 


= 2W + ((1 - 


. e MN)-n)/h B T\ N + 1 




+ s {^™ 


— ^mj 



(195) 



(196) 

Even though these equations are derived by neglecting all correlations, the stationary so- 
lution is given to the degree of approximation used — and large N — by firstly making the 
condensate chemical potential match that of the vapor, which requires to lowest order 

t*(N s ) = fi (197) 

and then requiring n^ = 0, which then gives 

1 



o£m/k B T J 
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and these are the solutions expected from statistical mechanics. The full time dependent 
solutions will be treated elsewhere, since they involve much more detail than we wish to give 
here. 



V. CONCLUSION 

A. The scope of this and future work 

We have provided in this paper a full quantum kinetic description of a Bose-Einstein 
condensate in contact with a bath consisting of a vapor of non-condensed atoms at a well 
defined temperature and chemical potential. The formulation up to Sect. Ill C, where the 
Bogoliubov approximation is introduced, is valid for any number of particles in the conden- 
sate, while the remainder of the paper, since it relies on the Bogoliubov approximation, is 
restricted to situations in which the number of particles in the condensate is rather large, 
perhaps of the order of a few hundreds. However an extension of the methodology to deal 
with small numbers of particles in the condensate would be straightforward, and requires us 
simply to express the field operators in terms of particle creation and destruction operators, 
and rewrite (50a-50e) appropriately. 

The key to the understanding of the condensate dynamics in this paper is the use of the 
particle number conserving Bogoliubov method [18], which enables us to use the Bogoliubov 
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spectrum as a function of the precise total number N of particles in the condensate band 
(and not merely the mean value (N)) , and further, to write precise expressions for the 
relevant transition matrix elements which appear in the master equation. The results of 
this procedure give very simple and easily understood equations of motion, especially when 
we; (a) ignore all quasiparticle effects; and (b) include only the simplest of the terms which 
contribute to the equations of motion, those which correspond to gain and loss from the 
condensate as a result of collisions between particles in the non-condensate band. This will 
clearly be possible at low enough temperatures and with sufficiently large condensates. The 
resulting equations are very like those used in laser theory [27-29], and thus are able to 
be handled by similar means. We have exhibited the basic theory of condensate number 
fluctuations; our treatment of the theory of phase fluctuations is simplified considerably, 
with a full treatment to appear in another paper [13], since these are measurable and of 
great interest, especially for the construction of an atom laser [38], and thus merit a very 
full and careful treatment. We will note here only that the methodology is straightforward, 
and very similar to that used for laser theory. 

The growth equation, (168), which represents in simplest terms the growth of the con- 
densate from a vapor has already been reported [26], and is in good general agreement with 
experiment. However, the process of condensate growth spans a wide range of computational 
and theoretical regimes, and the growth equation cannot be valid in all of these, since there 
will be corrections arising from both the inclusion of quasiparticle effects as well as those 
which must be made because the Bogoliubov approximation cannot be valid in the regime 
of small condensate numbers. The details of the modifications resulting from these effects, 
as well as those which arise because of depletion of the vapor, together with some detailed 
consideration of the experimental regimes will be reported in another paper [16]. 

The depletion of the vapor in the non-condensate band is an issue which requires for its 
confident handling a full treatment of the coupled kinetics of the condensate and the non- 
condensate band, and this will be reported in another paper [13], though the methodology 
will turn out mainly to depend on the use of equations of the Uehling-Uhlenbeck [15] kind 
for the non condensate, coupled to master equations, of the kind derived in this paper, for 
the condensate. 



B. Validity of the approximations 

The simplifications and approximations made in this paper can be summarized as follows: 

i) In going from (19) to (20) we have neglected the interaction between the condensate 
and the non-condensate bands in the kernel. This means that we assume that the 
condensate band has a spectrum which is not greatly changed by this interaction. 
For condensates presently being studied, this is certainly experimentally valid for the 
lowest levels, since their behaviour is well described by the time-dependent Gross- 
Pitaevskii equation. However, to give a quantitative criterion for the validity of this 
approximation seems at this stage to be too ambitious. 

ii) After the previous approximation, a Markov approximation is made. This amounts 
to requiring that the irreversible processes arising from the master equation are much 
slower than the reversible processes described by the separate evolution operators of the 
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condensate and non-condensate bands. For the non-condensate band, the arguments 
for its validity are the same as those in QKI, while for the condensate band, the 
evidence is experimental, that the damping measured is definitely weak compared to 
the timescale of oscillations, which can only be the consequence of the condensate 
band Hamiltonian Hq. 

iii) Another major approximation is the rotating wave or random phase approximation. 
Again, it is difficult to give a quantitative measure of its validity, other than to say 
that it should be valid when the basic behavior of the noncondensed vapor should be 
describable as transitions between states which are well described by single particle 
energy levels. All the exprimental evidence, including the very design of the experi- 
ments, seems to verify this. The basic premise of the random phase approximation is 
that the non-resonant terms are all of random signs, and effectively change sign fre- 
quently within times of interest, and it seems to us that the arguments for the validity 
of the random phase approximation are as good here as in similar situations in solid 
state physics. However, one can only make a definite statement by comparison with 
experiment. 

The question of validity can also be looked at from another point of view; namely, there will 
be traps, of the type illustrated in Fig. 1(a), in which the separation between the dynamics 
of the condensate and non-condensate band is as good as one wishes, and for which the 
approximations employed will certainly be valid. Even so, we believe that it is very likely 
that our approximations are already good for the Harmonic traps currently in use, but 
ultimately it seems to us that other kinds of trap may be both more useful, and more 
adapted to our kind of description. 

C. Relation to other work 

Apart from the work of Anglin [34], previous work on the field of condensate growth 
[30-33] has focussed on the growth of a condensate in a spatially homogeneous system, or 
an approximately spatially homogeneous system, rather than in a trapped system, and thus 
has not been applicable to the condensates presently being produced. The major difference 
in point of view in this paper is the description of the condensate in terms of the trap energy 
levels for a given number of particles, which is treated here by the modified Bogoliubov 
method, and the clear division into the condensate and non-condensate bands. It is these 
two innovations which make our treatment of condensate growth possible. 

Our simplified master equation (147) bears a very strong resemblance to that of Anglin 
[34] , who presented a model calculation based on a single level in a very tight potential well, 
centered inside a spherically symmetric bath of atoms moving essentially freely, apart from 
collisions. Anglin's model can be viewed as the description of the most ideal trap imaginble, 
in which the small trap permits only one energy level, and the larger trap is actually almost 
infinite. It is thus drastically simplified from the beginning, unlike ours, which treats all 
relevant features as exactly as possible. The term in Anglin's model which gives pure phase 
diffusion is not present in (147), but clearly arises from similar approximations applied to the 
term (142). Indeed if one takes the diagonal matrix elements of Anglin's master equation, 
one obtains an equation of the form (154) in which W + and W~ are related to each other 
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by (152), and the absolute value of W + is given by multiplying our value (153) (apart from 
the final term in curly brackets) by a factor of order of magnitude kT/Ef,, where Ej, is the 
binding energy of Anglin's potential as defined in his paper. 

The major assumption behind Anglin's derivation of the Markov approximation is that 
exp(/x//cT) is small. This requires that // be negative, and that — E\, < /i. Thus the physical 
situation is rather different from what we consider — the vapor of uncondensed atoms is 
always essentially classical, and the condensation only occurs because of the depth of the 
attractive potential. The collisions which establish the condensate are of the same kind as 
ours, but in Anglin's case a typical collision must yield one final atom in the condensate, 
whose negative energy ps —Ej, must be balanced by transferring the equivalent amount to the 
other atom, which must also take up all the momentum of the initial two vapor particles — 
this will involve typically particles with energy ps E^. If kT <C E\, this will be a rather small 
fraction of the total number of particles, and hence the extra factor kT jE^ seems quite 
understandable. 

Thus we conclude that Anglin's model is consistent with our analysis. 

D. The broken symmetry picture 

Our picture does appear to be at variance with the description of the condensate in 
terms of "broken gauge symmetry" or "Bose broken symmetry", which has been customary 
for many years [35-37]. It appears to us that these concepts are in fact not helpful in 
the description of the condensates which are produced in current experiments, and indeed 
we doubt whether they have any relevance to the physics. The stationary states of the 
system are indeed in our formulation simple statistical mixtures of number states of N, the 
total number of particles in the condensate band. The correlation functions (0'(x)0(x')) 
however do exhibit the mandatory property of "off diagonal long range order" , even though 
the stationary averages (0'(x)), (0(x')) are zero. This happens in laser theory as well, 
where the existence of a well defined phase of the optical field is the result of a rather long 
phase diffusion time compared with the other relaxation times — we will demonstrate the 
analogous property for the Bose condensate in [13], though the essential points have already 
been foreshadowed by several authors [39-42], who have shown how the measurement of 
the relative phase between two condensates introduces a relative coherence between the 
two condensates, which remains in existence on the time scale of phase diffusion. In a 
ferromagnet the concept of spontaneous symmetry breaking acquires its validity from the 
extremely long diffusion time for the magnetic moment of a sample of the ferromagnet — over 
a sufficiently long time the average magnetic moment would be zero, but this time is much 
longer than the lifetime of the average experimenter, who therefore does not see things this 
way. In the end it is a matter of taste where we draw the line between dephasing time which 
is very long and one which is essentially infinite — it depends on the relevant timescales of 
the experiment, including ultimately the timescale of the patience of the experimenter. 
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APPENDIX A: ESTIMATION OF THE SIZE OF THE CONDENSATE BAND 

In this appendix we will use the Bogoliubov method as described in Sect. Ill C to get an 
estimate of the energy Er above which we may neglect the effect of the condensate on the 
trap energy levels. 

a. Thomas-Fermi approximation 

We use the Thomas-Fermi approximation to the condensate wavefunction in which we 
neglect the Laplacian term in the time independent Gross-Pit aevskii equation, giving the 
approximate solution 



For a harmonic potential 

V T (x) = ax 2 + by 2 + cz 2 (A2) 

normalization requires 

2/5 

" (iV) = (— sr-j (A3) 

where 

a = mcul/2; b = muj 2 y /2] c = mu 2 z /2; (A4) 

b. Quasiparticle spectrum 
Use the definition (79) of 7i 3 , and substitute (Al) to get 

H, = /</ 3 xj - |l x t (x)V \-(x) + xHx)V TF (x)x(x)} 

+ f /^:(X'(X) + X(*)) ^ W ;^ W ) 

v(x)<Miv) V ' ' ' 

+ const. (A5) 
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with 

K TF (x) = max{V T (x) - fi(N), 0} (A6) 

We can then write 

H 3 = Y^ E q a\a q + J2 G kq {2a\a k + a\a\ + a q a k ) (A7) 

1 k,q 

where 

G kq = f d 3 x£,(x)a(x)UiV) _ y T (x)), (A8) 

E q are the energy levels of Vtf, and £ g (x) are the corresponding eigenfunctions, which have 
been chosen as real. 

c. Approximate evaluation of integrals 

For high lying excitations, we expect that the smooth nature of the ground state wave- 
function will mean that essentially G kq ~ S kq g q , so that the individual modes are diagonalized 
by a squeezing transformation 

a q = c q A q - s q A\ (A9) 

with, for small g q /E q , 

s q « g q /E q (A10) 

c q « 1 - S J (All) 

To estimate (? 9 , we approximate the wavefunction £ q by a WKB approximation to an ap- 
propriate three dimensional harmonic oscillator wavefunction (obtained by neglecting the 
flattening in the potential caused by the condensate, but not the chemical potential shift); 
this means we can approximate by smoothing over the oscillations, and then neglecting the 
x dependence, to get 

£ g (x) 2 « D ni ,n 2 ,n 3 = P x {n u 0)P y {n 2 , 0)P z {n 3 , 0) (A12) 

with 

1-1/2 



1 



P x (n u x) = - 
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(2m + i)h 

X 



muj 



X 



(A13) 



with similar equations for the other terms. This means that (letting q — »■ {711,712,713}, 

9n u n 2 ,n 3 = Aii,n 2 ,n 3 / rf 3 x(/i(A^) - Vt(x) 



'V(x)< M (AT) 

l6fx(N)V2 



/W 



15tt 2 \^ 3 n-=i(2n l + lV l 



(A14) 
(A15) 



n 3 

^n 1 ,n 2 ,n 3 ~^E( 2 ^ + 1 K-A t ( Ar ) ( A16 ) 
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d. Numerical comparisons 

For the case of the TOP trap of JILA using Rubidium, we find that the ratio 
9n 1 ,n 2 ,n 3 /E ni ^ n2tn:i can be chosen to be less than about 0.03 for all 711,712,713 on a fixed en- 
ergy surface provided E/fi(N) is more than about 2.4 at iV = 1,000 rising to 3.4 when 
N = 1,000,000. 

APPENDIX B: THE LOCAL EQUILIBRIUM APPROXIMATION FOR R NC 

We want to calculate 

u 2 



x5(K! + K 2 -K 3 -k) 

x5 (A 123 (0) - n{N)/%) F(K 1; 0)F(K 2 , 0) 

x[l + F(K 3 ,0)]|f(k)| 2 . (Bl) 

We assume 

i) That we can use the Maxwell-Boltzmann form 

F(K, 0) » exp r ^ j (B2) 

i) That -F(K3, 0) is negligible compared to 1 in the region of interest. 

ii) That the range of k in which |£(k)| 2 is non-negligible is very small. This means we 
can drop all the remaining k-dependence, and integrate |£(k)| 2 to give 1. 

We then carry out the K 3 integration to get 



W 



v 2 



(2tt) 5 ^ 2 

h 2 



d 3 Ki f d 3 Ki 



x5 



0m K K 2 1+ Kl- (K 1 + K 2 ) 2 ) - n(N) 



/n 



xexpp^M±Ml/^ (B3 ) 

From this point no further approximations need to be made. By changing variable to 
Q = K4 + K 2 , q = K4 — K 2 the integral becomes spherically symmetric in both of the 
variables, and the q integral can be carried out using the delta function, and the integral 
finally reduced to 

, u 2 2ir 2 m (Amn{N)\ 2 /2//-//(JV)\ 

w = n¥T" [-1^) exp [ k B T ) 

x/ ^" 1/2 (H-0 1/2 exp(-^), (B4) 
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Using [43], formula 13.2.5, the integral can be reduced to a form involving the confluent 
hypergeometric function C/(3/2, 3, z), which itself by formula 13.6.21 can be reduced to an 
expression involving the modified Bessel function Ki(z), which gives the final result 



Am {ak B T) 2 2ll/kgT J _/^y__ / fi N 
irh 3 \k B T l \k B T 



W+(N) = ^E^eW )™ Kl (J2L) , (B5) 



as in (153). 
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FIGURES 



(a) 





FIG. 1. The condensate and non-condensate bands for (a) an ideal trap and (b) a harmonic 
trap. 
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FIG. 2. Transitions in terms of quasiparticles 
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FIG. 3. Condensate growth as given by the growth equation for (a) and (b) Sodium. Scattering 
lengths are 5.71 nm and 2.75 nm, respectively. 
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